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A transport equation for eddy viscosity 

By P. A. Durbin 1 AND Z. Yang 2 


A transport equation for eddy viscosity is proposed for wall bounded turbulen 
flows The proposed model reduces to a quasi-homogeneous form far from surfaces. 
Near to a surface, the nonhomogeneous effect of the wall is modeled by an elliptic 
relaxation model. All the model terms are expressed in local variables and are coor- 
dinate independent; the model is intended to be used in complex flows. Turbulent 
channel flow and turbulent boundary layer flows with/without pressure gradient are 
calculated using the present model. Comparisons between model calculations and 
direct numerical simulation or experimental data show good agreement. 


1. Introduction 

Algebraic turbulence models (Baldwin-Lomax 1978, Cebeci-Smith 1974) have 
been used extensively in calculations of aerodynamic flows. These algebraic models 
are easy to implement numerically and give accurate predictions for simple flows, 
such as that over an airfoil with an attached boundary layer; however, they are 
inadequate when used for more complex flows, such as that over an air oi wi 
separation. In addition, these models contain nonlocal parameters in their formu- 
lations. The interpretation of these parameters becomes ambiguous in complex 
geometries. Higher order turbulence models have been proposed to overcome defi- 
ciencies of these algebraic models, the most popular among these being the k-t 
model. The empirical constants in the k - e model have been optimized in a variety 
of simple shear flows, giving what is commonly referred to as the Standard k e 

model (Launder and Spalding 1974, Rodi 1980). , . 

The Standard k - e model was devised for turbulence far away from boundaries 
For near wall turbulence, the model either has been used in conjunction with wall 
functions or has been modified by damping functions and other near wall terms to 
give the so called near wall version of the k - e model. This fixed up model can be 

integrated down to the wall. , . 

In the wall function approach, one assumes the existence of a universal wall layer, 
which is not valid in many complex flows. In the damping function approach the 
results depend on the damping functions used— a wide variety of which have been 
proposed. These damping functions themselves must be assumed to be of a universal 

It is also found that for near wall flows the k - e model is numerically stiff: it 
requires many grid points near the wall to get a grid independent solution. This 
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numerical stiffness makes the model rather unappealing for aerodynamic computa- 
tions. The numerical stiffness of the k e model is partly due to the rapid variation 
of e (and to a less extend k ) near the wall. 

The eddy viscosity, v t , varies much more gradually near the wall. This observation 
suggests that an equation for the eddy viscosity might be less stiff; this was one 
motivation for the work of Baldwin and Barth (1990) and Spalart and Allmaras 
(1992). In these papers, a transport equation for eddy viscosity of the type originally 
proposed by Nee and Kovasznay (1968) was formulated and solved in conjunction 
with the mean field equations. Wall effects were introduced in these models by 
damping functions, in which the distance to the wall y enters as an parameter. 
This wall distance y is not a local property; its definition may be ambiguous for 
flows in complex geometry — corner flow for example. Also, the damping functions 
are assumed to be universal, which does not seem valid in flows with strong adverse 
pressure gradients and separation. 

Modeling of near wall turbulence is very important for engineering calculation 
because the near wall portion of the boundary layer contributes substantially to 
the total momentum and heat transfer. In the near wall region, the turbulence is 
strongly inhomogeneous and strongly anisotropic. Recently, Durbin (1991) intro- 
duced the method of elliptic relaxation to model the blocking effect of the wall on the 
turbulence. In this method, the behavior of the near wall turbulence is given by the 
solution of an (elliptic) differential equation rather than by prescribed algebraic ex- 
pressions. Geometric information enters through the boundary conditions enforced 
at the wall. Thus, it is hoped that this formulation can handle flows in complex 
geometry. This approach also has the appeal that it is coordinate independent. 

In the present work, we propose a new transport model for the eddy viscosity. Far 
away from the wall, the model is based on the quasi-homogeneous approximation; 
near the surface, the wall effect is introduced via an elliptic relaxation equation. The 
structure of the paper is as follows: the proposed model is presented in section 2; in 
section 3, we show calculations of turbulent channel flow and turbulent boundary 
layer flows with and without pressure gradient; section 4 concludes the paper. 

2. The proposed model 

2.1 The eddy viscosity 

In the framework of the eddy viscosity model, the unknown Reynolds stress —HiuJ 
in the mean momentum equation is assumed to be related to the mean velocity field 
by 

2 

-UiUj =2u t Sij- -k6ij ( 1 ) 

where S{j is the strain rate of the mean field and v t is the eddy viscosity. The last 
term in (1), containing the turbulent kinetic energy k , can be absorbed into the 
mean pressure. Thus, only the eddy viscosity needs to be modeled. 

We propose that the general form of a transport equation for the eddy viscosity 
be 

(v + Vj/,1 + $ 


D 

Di Vt = V - 


( 2 ) 
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where the turbulent self transport is assumed to be analogous to the laminar dif- 
fusion and $ is a source term representing the combined effect of production an 
dissipation of u t . Our objective is to propose a model for * which is valid for both 
near wall and free shear layers. The justification for using a scalar eddy viscos- 
ity, as in (1), is that our attention is on thin shear layers, in which transport is 
predominantly transverse to the shear . 

2.2 Model for flow away from the wall 


Far from the solid surface, the inhomogeneity of a turbulence field is relatively 
weak. The model for $ could then be found by expanding about the homogeneous 

state. This quasi-homogeneous case is considered first. 

For homogeneous shear flow, if v t is much larger than v, the only variables that can 
enter the problem are the eddy viscosity v t and the shear rate 5. (The shear rate 5 
will be defined here via the rate of strain tensor of the mean field: S -2{b i} b a) ■) 
From dimensional reasoning, we have 


$ = $(S, v t ) - c\Sv t 


(3) 


where Cj is a model constant. Experimental results for homogeneous shear flow give 

Cl In general, the flow is not homogeneous. Flow inhomogeneity can be represented 
by IVi/fl and by |VS|, which are measures of the inhomogeneity of the turbulence 
field and of the mean field, respectively. If the inhomogeneity is weak, expansion 
about the homogeneous state suggests the form 

v 2 

$ = c\Svt — C2lVi/ t | 2 — C3-771 ( 4 ) 


where the length scale L„ is given by 



Other terms, such as \Vv t \v t /L v might be considered; (4) is analogous to the forms 
selected by Baldwin and Barth (1990) and Spalart and Almaras (1992). The model 
(4) allows inhomogeneities in both the turbulent field and the mean field to con- 
tribute to ‘dissipation’ of the eddy viscosity. 

2.8 Model for near wall turbulence 

Eq (4) is based on an expansion near the homogeneous state. The expansion 
is valid only for weak inhomogeneity. In the near wall region the flow changes 
very rapidly to adjust to the boundary condition at the wall; thus, the flow field 
is strongly inhomogeneous, and the expansion about the homogeneous state is not 
expected to be valid. (Direct application of equation (2) with $ given by equation 
(4) to turbulent channel flow at Re T = 395 leads to a value of the skin friction 
coefficient too high by about 30%.) In the present work, this strong inhomogeneity 
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is modeled by the elliptic relaxation model first proposed by Durbin (1991) for 
Reynolds stress closures. 

In the Reynolds stress transport equation, the velocity-pressure gradient correla- 
tion term is a good candidate for this type of relaxation because pressure is elliptic in 
nature. In the transport equation for the eddy viscosity, the pressure term does not 
appear explicitly; which term should be operated on by elliptic relaxation is ques- 
tionable. An examination of the transport equation for -uv/S in simple shear flows 
reveals that there are two origins for the elliptic relaxation: the velocity-pressure 
gradient correlation and the representation of v 2 in terms of uv. The combination 
of these two effects suggests, in the case of the eddy viscosity transport model, that 
we write cj = c u — c \2 and introduce the elliptic relaxation such that the transport 
equation appears as 



where P u is the quantity subject to relaxation. It is governed by the elliptic equation 
L 2 p V 2 P u - P v = -(c u Sv t - c 2 |V^| 2 ). (7) 

In (7), L p is the blocking length scale and is given by 

l\ = c\ maxQ,c, 2 ^). ( 8 ) 

The first term in the ‘max function represents the turbulent length scale, which 
is appropriate away from the wall; the second term is chosen very close to the wall, 
where the appropriate length is the Kolmogorov scale. The switching from one 
length scale to the other is controlled by the value of c/ and the total length scale 
is controlled by c p . 

2.4 Modification near 5 = 0 

Equations (5) and (8) give expressions for the length scales. These estimates 
break down when 5 = 0. 5 = 0 can occur at the edge of a boundary layer or inside 
a turbulent field, for example at the centerline of channel flow. A revision to the 
length scale formulation is needed to accommodate these cases. The singularity will 
be removed by the following modification: 



Over most of the flow field, the first term on the right hand sides of (9) and (10) 
is much larger than the second term. The modification only affects the region near 

5 = 0. 
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In the course of developing and testing the model during the summer program, 
another modification to L v was inadvertently introduced: in the outermost portion 
of boundary layers, L v became constant. The algorithm was to make L„ be in- 
dependent of y when dv t /dy < 0. u t decreases when y/6 99 is greater than ~ 0.9. 
In future work, this restriction will be removed; it is an effective way to avoid the 
difficulty of specifying L v in the free-stream, but its applicability is restricted to 
attached boundary layers. 

2.5 Boundary conditions 

Equations (6) and (7), together with the length scales given by (9) and (10), 
are the eddy viscosity model we are proposing for wall bounded turbulent flows. 
They form a fourth order system of equations, and, consequently, four boundary 
conditions are needed. These conditions are as follows: On a solid surface 

v t = n • = 0, (H) 

where n is the surface normal; at a symmetry plane, the center line of a channel for 

exampk ’ „; = n = o, <>2) 

where the prime denotes the derivative normal to the symmetry plane; at the 
freestream edge of a boundary layer or shear layer 

v' t = P v = 0. ( 13 ) 

The prime in (13) represents differentiation normal to the edge of the boundary 
layer. In some situations, (13) might be replaced by a condition prescribing prop- 
erties of free-stream turbulence if it is present. 

2.6 Model constants 

There are seven empirical constants in the present model. They are a, c lu c n, 
c 2 , c 3 , c p , ci. c u and c J2 are related by 


cn — CJ2 = C\. 


ci ss 0.12 was determined by experiments on homogeneous shear flow. 

Another relation among the constants is found by requiring that the model give 
the right behavior in the logarithmic region of a boundary layer. This imposes the 
constraint 


where /c is the Von Karman constant ; we take k 0.41. 

The constant a in the turbulent diffusion term should be about 1 and the value 
of Cp, which marks the switching from the Kolmogorov length scale to the turbulent 
length scale, should be about 3.5. All the other constants are chosen by comparing 
with direct numerical simulation data for turbulent channel flow at Re T = 395 and 
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with experimental data for a zero pressure gradient, flat plate boundary. The values 
of the constants selected in this study are 

<* = L3 > c n = 0.4, c \2 = 0.28, c 2 = 1.3, c 3 = 0.2, c p = 3.5, c t = 1.6. (14) 

The constants listed above represent a preliminary choice: they are subject to 
adjustment in the course of further model development. 

3. Results and discussions 

Turbulent channel flows at different Reynolds numbers and turbulent boundary 
layers with and without pressure gradient were calculated using the present model. 
An implicit finite difference scheme was used to solve the momentum equation and 
the turbulence equations. For the present cases, the equations are parabolic and 
a marching scheme was used. The finite difference equations for v t and P v were 
solved as a coupled, block tridiagonal system. Variable grid spacing was used, with 
densest spacing near the wall. The total number of grid points was set to 105, which 
is large enough for the results to be grid independent. 

Two dimensional, fully developed channel flow is attractive for model testing be- 
cause it is statistically steady and nonhomogeneous in only one direction. Solutions 
can be found very efficiently, yet the effects of the wail on turbulent shear flow 
are still present. Computations were carried out for 2D fully developed turbulent 
channel flows at Re r = 180 and Re r = 395, for which DNS data are available 
(J. Kim, private communication). Figure 1 shows profiles of the mean velocity at 
these two Reynolds numbers along with DNS data. Both the dependent variable 
and the independent variables are represented in wall units. The predictions are in 
reasonable agreement with the data. 

The second example calculated is the turbulent boundary layer with zero pressure 
gradient. The zero pressure gradient boundary layer on a flat plate has a self-similar 
profile (in the outer portion of the boundary layer, when scaled on outer variables). 
Thus, arbitrary profiles could be used as the initial conditions, and the solution 
should develop into its similarity form. Figure 2 shows velocity profiles scaled 
by the free stream velocity and the boundary layer thickness at Re $ =3195 and 
Ree = 5473, respectively. The experimental data is from Coles and Hirst (1968). 
It is found that indeed the velocity develops into a similarity profile, which serves 
as a good check for the present model. Figure 3 contains a conventional log-linear 
plot of the mean velocity profile. The agreement with the data is excellent except 
in the viscous sublayer, where the model profile is a bit low. 

From an engineering point of view, the skin friction is of great interest. The skin 
friction coefficient as a function of Re $ is shown in figure 4, along with experimental 
data (Coles and Hirst 1968). It is seen that the model gives a good prediction 
of C /, The growth of the boundary layer in this range of Reynolds number is 
illustrated by figure 5, which shows displacement thickness as a function of Re*. 6* 
is normalized by the boundary layer thickness at Re$ = 1000. The model gives a 
slight overprediction of the boundary layer growth rate. 
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FIGURE 1. Mean velocity profiles in channel flow. •, Re r — 395; x, Re T 180. 
The lower Reynolds number profile is offset 5 units for clarity. 


Next we calculate the development of a turbulent boundary layer in a pressure 
gradient. The test case chosen is the Samuel and Joubert (1974) experiment on 
a boundary layer developing into an increasingly adverse pressure gradient. This 
boundary layer is not self-similar. During the 1981-82 Stanford Conference on 
Complex Turbulent Flows, it was found that this flow is very difficult to predict; 
since then, it has become a standard test case for turbulence models. 

The initial condition for the computation was set in the following manner: We 
assume that the turbulent boundary layer develops under zero pressure gradient 
up to a point, x, say, at which the pressure gradient (which is known from the 
experiment) is applied. The position i, is determined so that at the first point 
of the working section of wind tunnel, where the experimental data begin, the 
predicted values of Re e and Cf agree with those of the experiment. 

The predicted development of the boundary layer is shown in figure 6 along with 
the experimental data. The variation of the skin friction coefficient with x and 
the growth of the boundary layer thickness are predicted very well. The velocity 
profiles at x = 1.87 m and x = 2.55 m are shown in figure 7. The agreement with 
the experimental data is again found to be reasonable. 
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Figure 2. Similarity form of mean velocity profiles in zero pressure gradient 
boundary layer. •( ), Reg = 3195 ; x( ), Reg = 5473 . 



Figure 3. Log-linear plot of mean velocity profile in zero pressure gradient 
boundary layer. 
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FIGURE 4. Skin friction coefficient versus momentum thickness Reynolds number 
for zero pressure gradient boundary layer. 



FIGURE 5. Boundary layer thickness versus Reynolds number for zero pressure 
gradient boundary layer. 

4. Conclusions 

We have presented a transport equation for eddy viscosity for wall bounded flows. 
Away from the wall, the proposed model reduces to a quasi-homogeneous model. 
The near wall effect is represented by an elliptic relaxation equation. Since geomet- 
rical information comes through the location where the boundary conditions are 
enforced while the model itself is free from the boundary information, the proposed 
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FIGURE 6. Skin friction coefficient and displacement thickness in the Samuel and 
Joubert boundary layer. 



Figure 7. Mean velocity profiles at x - 1.87m and 2.55 m in the Samuel and 
Joubert boundary layer. 
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model has the potential to handle flows in complex geometry. All the quantities 

used in the model are local properties and are coordinate independent; the model 

is also Galilean invariant. . 

Channel flow at two Reynolds numbers and flat plate turbulent boundary layers 

with zero and adverse pressure gradients were calculated using the present model. 

The comparison with the DNS and experimental data were found to be quite promis- 
ing. 
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IV. The combustion group 


The combustion group conducted six projects. Three projects were related to 
premixed and three to non-premixed combustion. The invited participants were: 
Mr. M. Baum (Ecole Centrale Paris), Dr. J. H. Chen (Sandia National Labora- 
tories); Prof. R. 0. Fox (Kansas State University), Dr. D. C. Haworth (General 
Motors Research Laboratories), Prof. J. C. Hill (Iowa State University), Prof. S. 
Mahalingam (University of Colorado), Dr. T. Poinsot (Institut de Mecanique des 
Fluides de Toulouse), Prof. I. K. Puri (University of Illinois). The local partici- 
pants were: Dr. R. D. Moser and Dr. M. M. Rogers from NASA Ames Research 
Center, and Dr. F. Gao, Dr. A. Trouve, and Dr. L. Vervisch from the Center for 
Turbulence Research. 

The broad scientific objectives of this group were similar to those of the CTR 
1990 Summer Program: understanding of fundamental phenomena controlling tur- 
bulent combustion and application to modeling. However, the tools used in 1992 
have covered a wider range: 2D or 3D, variable or constant density, simple or com- 
plex chemistry formulations for Direct Numerical Simulation (DNS) of turbulent 
combustion have been used. Recent progress in numerical analysis and code de- 
velopment have allowed us to use tools which were well adapted to the physical 
problems considered by each group. 

The first three projects were aimed at increasing our understanding of turbulent 
premixed flames. 

Poinsot k Haworth studied the interaction between a turbulent flame and a cold 
wall. Flame quenching distances as well as wall heat fluxes were measured from 
2D simulations. The characteristics of flamelets reaching the wall (curvature, flame 
speed, quenching times) were used to build a ‘law-of-the-wall’ model for turbulent 
combustion. A simplified version of this model was derived and may be implemented 
in any flamelet model for turbulent premixed combustion. 

A new 3D variable-density code was used by Trouve k Poinsot to investigate 
the modeling of the evolution equation for the flame surface density. This analysis 
shows the limits of present flamelet models and suggests how the exact evolution 
equation for reactive surfaces may be closed to provide a suitable model. The effect 
of the Lewis number was evidenced and its influence on source or consumption 
terms for the flame surface density was demonstrated. 

The first objective of the work by Baum, et al. was to prove the feasibility 
of DNS with complex chemistry and the potential of this approach for pollution 
studies. A second goal was to check whether the DNS results previously obtained 
with single-step chemistry are modified by accounting for more realistic chemical 
schemes. The project was based on a new 2D code where a DNS technique was 
coupled to CHEMKIN and TRANSPORT, the SANDIA packages for reacting flows 
with complex chemistry. Using this tool, it was possible to investigate the structure 
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of H2-02 turbulent flames with the Warnatz scheme (9 species - 19 reactions). This 
project also provides the first analysis of the effects of stretch and curvature on 
flames with realistic chemistry. 

The first project of Chen et al. was related to the effects of finite rate chemistry 
and differential diffusion on the structure of turbulent non-premixed flames. 3D 
variable- density simulations were performed over a range of conditions correspond- 
ing to feist and slow chemistry. This project provided new insights on the validity 
of the flamelet assumption as well as on the effects of transient regimes and small 
scales on the inner structure of the flame zone. 

Chen et al. investigated one of the classical assumptions used to model turbulent 
non-premixed flames, i.e. single-step chemistry. Single-step and two-step chemical 
schemes were compared using 2D variable-density simulations of a non-premixed 
flame in isotropic turbulent flow. An important result was that extinction lim- 
its appear to be quite different: while single-step chemistry lead to multiple local 
extinctions, two-step chemistry feature more robust flames which do not quench. 

DNS of a single-step chemical reaction with non-premixed reactants in forced 
isotropic turbulence were used by Fox et al. to obtain joint pdf’s and other statistical 
information to parameterize and test a Fokker-Planck turbulent mixing model. The 
simulations were performed using a constant density, 3D spectral code developed 
by Moser and Rogers. Physical features as well as various statistics of the reacting 
scalars and their gradients were examined and compared to the model. 

This Summer Program has brought many new and original results. The activity 
on non-premixed combustion was more intense than during previous programs and 
opened new perspectives for modeling for those flames. The demonstration that 
DNS of reacting flows was possible while taking into account complex chemistry or 
the presence of walls also opens new fields of investigations. 

Modeling was one of our first objectives in this work. A fundamental aspect of the 
1992 work is the impact and the power of DNS to answer certain critical questions 
for turbulent combustion models. We believe that DNS of reacting flows is now 
reaching a point where individual terms in combustion models may be estimated 
from DNS (as done by Trouve &: Poinsot) and that this possibility will change the 
way we construct models in the next few years. 


Thierry J. Poinsot 
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